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ABSTRACT 

A numerical lunar theory which can be used to obtain the 
rectangular coordinates for a satellite moving in a highly inclined 
orbital plane is developed. The arguments of the theory are the 
linear functions of the true orbital longitude of the satellite and of 
the sun and thus are of Laplacian type. Hansen's device is em- 
ployed to perform the integration and for this purpose a fictitious 
satellite is introduced whose true orbital longitude is considered 
as a constant until the integration is completed. The perturba- 
tions in the orbital plane are obtained by means of a W- function 
analogous to the W-funetion of the classical Hansen theory and 
the perturbations of the orbital plane carried by the vector r , the 
unit vector of the fictitious satellite relative to the mean orbital 
plane. It is shown here that the combination of ideas of Laplace, 
Hansen, and Hill represents a convenient way to obtain a numeri- 
cal lunar theory. At the same time, this work also represents a 
further development and a simplification of the results given by 
the author in his previous article. 
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ON THE NUMERICAL THEORY OF SATELLITES 
WITH HIGHLY INCLINED ORBITS 


by 

Peter Musen 

Goddard Space Flight Center 


INTRODUCTION 

In this article, we develop a form of the differential equations of the lunar theory which can 
be used to obtain the rectangular coordinates for a satellite moving in an highly inclined orbital 
plane. The integration of these equations is based on solving a certain linear partial differential 
equation by means of iteration. The author suggested this idea in an earlier work (Reference 1). 
The arguments of the theory are the linear functions of the true orbital longitude of the satellite, 
v, and the true orbital longitude, v 4 , of the sun and thus are of Laplacian type. The first sugges- 
tion to use v in the theory satellites comes from Brown (Reference 2). The use of the true orbital 
longitude speeds up the convergence of the development of the disturbing function, compared to the 
development in terms of the mean anomalies. 

As in Hansen ! s theory (Reference 3), we split the perturbations of the satellite into the per- 
turbations in the orbital plane and the perturbations of the orbital plane. We make use of 
Hansen T s device to perform the integration and for this purpose introduce a fictitious satellite 
whose true orbital longitude, w, is considered as a temporary constant until the integration is com- 
pleted. After the completion of integration, we apply the "bar operation" by replacing w with v, the 
true orbital longitude of the real satellite. 

It is of interest to note that the W-function determining the perturbations in the orbit plane is 
simpler in the theory of Laplacian type than in the classical Hansen's theory. Thus the combina- 
tion of ideas of Laplace, Hansen, and Hill represents a convenient way to obtain a numerical theory 
of the satellite. The exposition given here represents also a further development and a simplifi- 
cation of the results given by the author in his previous article (Reference 4). In the classical 
Hansen theory use is made of three auxiliary parameters, 

P = 2 sin ^ i sin N' , 

Q - 2 sin ^ i cos N' , 


1 



and 


K' , 


where N' is a purely periodic part in -(0 +a)/ 2 and K' is a purely periodic part in +(0 ~a)/ 2 . For 
the sake of the symmetry the author (Reference 5) has suggested the use of four parameters, 


k, = sin — i cos N', k- = cos — i 
1 2 3 2 


K' 


sin — i sin N', 
2 


and 


A. 4 = cos — i cos K' 


$ 


The implicit role of Hansen’s, as well as the author’s parameters is only an auxiliary one: they 
only help to form the components of the unit vector of the satellite with respect to the mean orbit 
plane. For this reason we discarded their use and resorted to the determination of this unit vector 
directly. 


The differential equation governing the perturbations of the unit vector is simple enough to 
justify this modification and the total number of the differential equations is reduced for one. The 
method given here is equally applicable to the planetary satellites disturbed by the sun or to the 
lunar orbiter disturbed by the earth, as well as to the artificial satellites of the earth whose motion 
is disturbed by the presence of the zonal and the tesseral harmonics. In the first two instances, 
the eccentricity must be small or moderate and in the last instance, it does not include the case of 
critical inclination. 


BASIC DIFFERENTIAL EQUATIONS 

We refer the motion of the satellite to a moving ideal system of coordinates whose x and y 
axes lie in the osculating- orbit plane. The equations of motion can be written in the following 
standard form: 


d 2 r r /d v\ 2 1 

d7 2 " r \dt/ ’ ~ 7 2 + ^ (i) 


A ( r 2 _ AL ( 2 ) 

d t \ d t / 3 V 


Setting 


2 dv J_ 

d t h 


( 3 ) 
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and 


T ’ (4) 


we can write 


u - h 2 + h 2 e cos (v - x) . (5) 

The angle / can be decomposed into the secular part, 


(X) = (1 - gj) V + 7 T 0 


(6) 


and into a purely periodic part, <p. Thus, 


X - (1 - gj) v + 77 q + 4 > . 


(7) 


From Equations 5 and 7 we have 


u = h 2 + eh 2 COS (Vj - <fi) 


(8) 


where we set 


V 


1 


Si v - ^0 


(9) 


The mean value of u can be defined as 


u = h 0 + e o h o cos V v 

where e 0 is the constant of the eccentricity and l/h 0 is the constant of the n area integral Both ele- 
ments, e 0 and h Q , together with other elements, must be chosen in such a way that no secular or 
mixed terms appear in the development of the coordinates into trigonometric series. 

In order to set Hansen's integration procedure it is convenient in the Laplace- Hansen type of 
lunar theory to define the basic W- function by means of the equation, 

h h 0 

W = — [1 + e cos (Wj - 0)] (1 + e 0 cos , (H) 

h o " 

where w t is considered as a temporary constant. 
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In our case the Hansen bar operation consists of the replacement of w x by v x after the inte- 
gration is completed. In forming the differential equation for w, as well as in the process of integra- 
tion, w x is invariable. This form of the W- function is different from the classical one but it leads 
to a simpler differential equation for its determination. The application of the bar operator to 
Equation 11 gives, 


u = u + hj- W . 
h 0 

When the stretching factor, 1 + v, is introduced by means of the equation 


(12) 


u = (1 + y) u , 


we obtain from Equation 12 


1 + V = 



(12a) 


(13) 


The practical way of computing v by means of iteration is based on the use of the formula 


V 


(1 + v) £ • 
hn 




W. 


(13a) 


By use of the equation (Reference 3) 


_d_ 
d t 


-t- e cos (x - P) 


i an 

= F 0 3r 


s i n ( v - P) 


+ 


Bv 




COS ( V - p) 



e ^7 sin (x - P) , 

d t 


(14) 


and setting 


P = v - v x + w x 
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we obtain 


_d 
d t 


— e cos (Wj - 0) 
_ h o 


1 3fl . , . 

= r~ — sin (v| - w x ) 
h Q dr 


+ h o (— + — -) cos ( y i - w i^ 

3v \K K) 


+ (1 - g,) e ^7 sin ( 4 > - Wj) 

h o dt 


( 15 ) 


Combining Equation 15 with both of the standard equations 


d h h 2 an 

' n _ 


dt h 


o h 2 3v 


and 


we deduce that 


_d _ h — 
dt h + 0 Bv 


dW 

dt 


h 2 


L\ h 0 h 0 


+ \ C os (Vj - XV,) - — + 1 + e 0 cos 


h 2 




an 

Bv 


+ — sin (v. - w t ) + £- e (1 - gj) ^ sin (<£ - Wj) 
h Q 3r 1 1 h 0 1 dt 


(16) 


When Equation 12 is taken into account, we obtain 


dW 
d t 


h 2 u h v 

+ +— w cos (v. - w) 

h 2 h 2 \ ) 


h 2 - 

— - + 1 + e n cos w 


an 

Bv 


1 BQ . N h \ dv . . 

+ h~ 3^7 sin <^ v i “ w i> + — e C 1 - — sin (0 - Wi) 


dt 


(17) 
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When sin (0 - Wj) is eliminated from the last equation by means of the relation 


h BW 

— e sm (<f> - Wj) = — 
h n 1 o w , 



sin Wj 


( 18 ) 


and 


dt _ h 
dv ” u 2 


is taken into consideration, we have 


dw M bo „ h o Bn 

— = N — |M~ 

dv Bu u 2 Bv 


+ (1 - 


Bff 

Bwj 



sin w 


(19) 


where 


N = - A sin (v - w ) 
h o 

and 



'h 2 h — 

— + 1 + e n cos w t + — W 

X h o 


COS (Vj - Wj) 


+ 1 + e Q cos Wj 


The motion of the pericenter, V - g p is obtained from the condition that no term of the form A sin Wj 
appears in Equation 19. 

If the eccentricity, e 0 , is not very small, say approximately 0.01-0.3, then the numerical values 
of all the elements can be substituted from the outset. If the eccentricity is smaller than approxi- 
mately 0.01, then in order to avoid the numerical difficulty in the determination of the motion 
of the perigee, it is recommended that the first power of the eccentricity be kept in the literal 
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form. In other words, every trigonometric series, T, must be written in the form 

T = A 0 + e 0 Aj , 

where A 0 is independent from e 0 and both A 0 and A x are trigonometric series with purely numeri- 
cal coefficients. We write the product of two such series in the form 

(A q + e Q Aj) (B 0 + e Q Bj) = A q B q + e Q (A J B Q + A Q Bj + e Q Aj B x ) 

and the numerical value of e 0 is then substituted in the parentheses. Such a simple device serves 
as a safeguard against the small divisor in the determination of the motion of the pericenter from 
Equation 19. The unit vector r° can be written in the form 


7 ° = A a (0) • Vi) • A 3 (-o0 • 


COS V 

sin v 
0 


( 20 ) 


where 


and 


Aj(a) = 


+ 10 0 
0 + cos a - sin a 

0 +sina +cosa 


A 3 (a-) 


+ cos a - s in a 0 
+ s in a + cos a 0 

0 0+1 


Thus we have 


7° = A 3 (0) - Aj(i) • 


COS (v - or) 
sin (v - o) 
0 
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we can then set 


V - a = v 2 + N' + K' 


and 


e = v 3 - N' + k' , 

where N' and K' are purely periodic and v 2 and v 3 are the linear arguments with respect to v, such 
that 


and 


v 2 = V + ^0 


V 3 = v + 


(21) 


The argument v 2 represents the mean argument of the latitude and v 3 represents the mean longitude 
of the ascending node. The coefficient g 2 is of order of one and g 3 is of the order of perturbations. 

Using the auxiliary parameters 


1 . 


k, = sin — i cos N', 
1 9 


k? = cos - i sin K' , 
3 2 


k~ = sin — i sin N', 
2 2 


k A = cos -i i cos K' , 
4 2 


the unit vector 7° can be represented as (Reference 4) 


COS V 2 


r° - A 3 (v 3 ) * A 


sin v 2 


0 


(22) 


where the matrix A carries all the periodic effects in i and e. The elements \ of the matrix A 
are simple polynomials in k 19 k 29 k s , \ as: 


\i = +A?-A?-A|+^. 


X 21 = + 2 ^3 \ ~ 2 \ ^2 ' > 


^•31 — ^ k 3 A.j + 2 k 4 


(23) 
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\_12 - - 2 “ 2 \ 2 


^22 “ \ + ^2 ^3 + ^4 


X- 32 _ + 2 X.j X- 4 - 2 \ 2 \ 3 


\ 13 - + 2 \ 3 - 2 k 2 X . 4 


\ 23 — 2 \j \ 4 - 2 \ 2 \ 3 


v 33 - - - 2 r 3 T 4 


The author has established (Reference 6) the following set of the differential equations for the 
variations of the parameters: 


^ a - g 2 + g 3 ^^ 2 +^ h + (A i + x 3> - (X 1 + *a V 5 XJ ' (26) 


d\ 1 

— 1 = + i (1 - g 2 + g 3 

dt 2 2 3 


) £ + i h - <*3 + *1> ^ ^ - \ V ^ + (\ \, + ^ V ^1 > (27) 


d X-t i 

- 3 = + { a - g 2 - g 3 


^ ^ + ? h [" ^ + ^ ^ + ^ + ^ ^ + ^ ^ S; 




^ h j^+ (\j + + (>>. 2 \ ^ 3 ) - (\ \* + ^2 ^ 3 ^ 


From Equation 22-28 we deduce that 


[ + a-B2^ 2 + «3Si] £+i hA ^ 
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and 


^ 12 - [-a - g 2 ) + g 3 *- 22 ] 1 h B \ 3 


d t 


(31) 


^21 r XX x 1 dv 1 


d t 


[.(l-g 2 )X 22 -g 3 \ 1 ] 3T + o hA ^ 


(32) 


^22 r ,* xx x i dv 1 


dt 


[~ (1 ^21 “ g 3 ^12] + 2 ^ ^ 


(33) 


^31 
d t 


+ a-8 2 H 32 ^ + ihA\ 3 3 


(34) 


d\_ 


d t 


= - (1 - e 2 >^i p t ^ h B ^3 


(35) 


where we set 


A = - 


30 

3X7 


+ ^x 


an 


-K 


m . an 
bT +x » 


3\. 


(36) 


and 


30 


30 


30 


30 


+ ^3 + \l ^x ” N S ^ 


3\, 


3A. 


3\, 


(37) 


The simplicity of Equations 30-35, as well as Equation 22, show that we can discard the use of 
^-parameters and introduce the unit vector, 


r = Tj i + r 2 j + r 3 k, 

instead, where 

Pi = \ u cos W 2 + \ 12 sin W 2 > 
f 2 = \ 21 cos w 2 + \ 22 sin w 2 , 


(38) 


(39) 

(40) 
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(41) 


r 3 = \ 31 cos W 2 + k 32 sin W 2 , 


and 


r i + r a+r> = l . 

We shall consider w 2 constant during the formation of the differential equation for P and during the 
integration. After the integration is completed, we apply the bar operation and replace w 2 with v 2 . 
The vector p is the unit vector of the satellite relative to the mean orbital plane. Thus, 

f = T 1 I + r 2 T + r 3 k . (38a) 


where 

r x = \ n cos v 2 + \ 12 sin v 2 , (39a) 

P 2 = k 21 cos v 2 + k 22 sin v 2> (40a) 

r 3 = \ 31 cos v 2 + \ 32 sin v 2 , (41a) 

and 




r 2 + r ? = !■ 


and Equation 22 takes a concise form, 


r = a 3 (v 3 ) • r . (42) 

In performing the iteration we have to distinguish between v 2 in the development of perturba- 
tions and the elliptic v 2 in the vector r. The derivatives of the disturbing function are formed 
with respect to the elliptic v 2 . For this reason it is convenient to separate the elliptic and the 
non-elliptic v 2 in the disturbing function by temporarily replacing the vector r with P. The dis- 
turbing function, so modified, will be designated by n. Taking Equation 39 through 41 into account, 
we obtain 


1 / an 'dn .an , an \ ( an ^n_ x bi i_ x \ 

2 (-^3^ 13+ 3r 2 ^ 3 3 r 3 3 3 j 


sin w, = - R 


■ ^ II sin w 2 . 
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i / . an .an .an . an \ /an 

i = * IsfT 


. an . 

^13 + 3P - ^3 + 


an x \ 


cos w 2 = + R • Vp n cos w 2 , 


and 


- A = - C sin v„ 


and 


-B = + C cos v, 
2 


where 


r _ Bn , Bn v Bn , 

= Bf; Xl3 + Br 2 ^ 23 + Br 3 33 


The bar operation here means replacing w 2 with v 2 . Equations 30 through 35 become: 


dk. . -iHv_ 

[+ (1 - g 2 ) \ 2 + g 3 ^21 ] ^~ hCK l3 Sln V 2 


dt 

dk. 


12 


dv 


dt - [- (1 - g 2 ) + g 3 X 22 ] JjT + h CX 1 


d k 0 . n H V 

= [+ (1 - g 2 ^ ^22 - 83 TT- hC ^23 sin V 2 . 


dk. 

d t 


*2' 22 - & 3 'Ml J dt 

dv 


22 = [- (1 - g 2 ) k 21 - g 3 \ 2 ] ^ + h C \ 23 cos v 2 , 


d ^31 /, _ X . dv 


d t 


- + (1 - g 2 ) ^32 ^ - h C \ 33 sin v 2 , 


and 


^32 N v dv , * . 

= - (1 ~ i,) K, TT + h C ^33 COS V 2 


dt 


d t 


From the above equations we deduce that 


dr, Br h 2 h 2 

-r 1 = + (1 - g 2 ) aT- + g 3 r 2 +-T c *13 sin ( W 2 - V * 
dv dw 2 h 2 u 2 


(43) 


12 


(44) 


an 


o dl 9 h 2 

— = + (1 - g 2 ) - — - g 3 r + — ' — CL sin (w 2 - v ) 

d v ^ w 2 h 2 u 2 2 2 


and 


dr. 


3 s ^3 h 2 


- + (1 - g 2 ) - — + C \ sin (w 2 - V ) 


dv 


3w 2 h 2 U 2 


(45) 


or, in the vectorial form, 


dr 

dv 


Br 


h 2 h o 


h 2 - 


+ (1 - g 2 ) ^ — + g 3 r X k + — * — R C sin (w 2 - v 2 ) 
3w 2 h 2 u 2 


(46) 


The components \ 13 , k 23 , k 33 of the vector R are computed using the formula 


R = ( r * 


jr 

Bw, 


2 1 w 2 =0 


For the lunar case (taking Equation 42 into account), we obtain 

cos H = r ° • = [cos v 4 , sin v 4 , 0] ■ A 3 (v 3 ) • V 


= Tj cos (v 3 - v 4 ) - r 2 sin (v 3 - v 4 ) 


(47) 


Thus, in performing the transition from n to n we must substitute the expression 


Tj cos (v 3 - v 4 ) - r 2 sin (v 3 - v 4 ) 


for cos H. From the series 


— P, (cos H) + m ' 1-3 P, (cos H) 

r' 3 2 #4 3 


+ P 4 (COS H) + 


we obtain 


n = k 2 h 2 (i + v) 2 


/j?v ('l 
[h r2 J \ U 


3 p2 3 p2 _ 1 

4 1 4 2 2 


13 



I I III INI 


■I II III I 


II 


+ 1 ( r i - r 2 2 ) cos (2 V 3 - 2 v 4 ) - § r, r 2 sin (2 V 3 - 2 v 4 ) 


+ A. 2 h 2 a (1 + v) 


U 2\3 
3 



( + T r i + T r i r 2 “ | r i) COS (V 3 - v 4^ 


+ (-T r ?r a -£r| + §r a ) sin ( v 3 -v 4 ) 
+ ( + | r i "T ri r ^) COS ( 3v 3 ~ 3 v 4 ) 

+ (-¥ r i^ + f^) Si "< 3v 3- 3 V 


+ \ 2 h 2 a 2 (1 + i/) 4 



r 4 , ^ r 2 r 2 
1 + ~32 1 2 


105 p4 15p 2 15p 2 3\ 

2 T 1 8 2 + sj 


+ (+ — r? - — r 2 l 3 r? + rA cos (2v,.2vj 

^ 16 1 16 2 8 1 8 2 / 3 4 

+ (-f ri r 2-T r i r3 2 + T r ' r ?) -n(2v 3 _2v 4 ) 

4 ( + H r i + if r 2 ~ ^ r i 2 r 2 2 ) cos(4v 3 -4v 4 ) 


+ 





sin (4 v 3 - 4 v 4 ) 


+ 


14 



where 



k is the analogue of the parameter m of Delaunay 1 s theory, and a is the parallactic factor of our 
theory. Equation 19 becomes 


dw = N in + M ^o an 

dv 3u u 2 Bw 2 


+ (1 - g 2 ) 





sin w 


(19a) 


Equations 43 through 45 and Equation 19a can be simplified in the case of the artificial satel- 
lite, provided only zonal harmonics are considered. We have 


n = k 2 u 3 (i _ 3 r 2 ) + k 3 u 4 (3 r 3 _ 5 r 2 ) 


+ k 4 u 5 (3 - 30 r 2 3 + 35 r 4 ) + 


dr i ^ r i 

^= + a-g 2 )^ + g 3 r 2 


h 2 h o an . , v 

-f ^"33 ^13 ® i ^ (w 2 ~~ ^2^ 9 

h 2 U 2 3r 3 3 13 2 2. 


dr 2 3r 2 


+ — •—! \ \ sin (w, - vO, 

h 2 u 2 3r 3 


dr 3 3r 3 

d7= + (1 - g2) ^ 




h 2 


h o an >• 
75 3TT 


sin (w 2 - v 2 ), 


15 


and 


aw _ N Bn M ^o Bn_ . 

dv Bu + u 2 3r 3 Bw 2 


+ (1 - g 2 ) 




We start the iteration with 


h. = i, i + v = i, 
h o 

r = cos W 2 , T 2 = cos i Q sin w 2 , and T 3 = sin i Q sin w 2 

and repeat it until we reach the final values. The perturbed coordinates are the trigonometric 
series in four arguments: 


Vj, V 2 , V 3 , and V 4 . 

In the classical Laplacian theory the angle v 4 is eliminated in favor of the angle mv-c 4 , where c 4 
is a constant. We can avoid this elimination, as well as the substitution of trigonometric series 
into the arguments, if we base the solution on integration of a partial differential equation instead 
of performing the quadratures. 

We have 


^4 _ h_ M 2 
dv h' \u/ 


/l + m' , 


or 




V) 2 


h / u'/h' 


h oV u/hg 


(48) 


Designating the constant part in the right side of the last equation by g 4 , we have 


a a o o o v ° 

= g » * 7 , + 3^: + 63 + 64 3^: + K 3^7’ 
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where we set 


= (£Y.* a+ vW^V- 

\ h o/ h o \ u/Uq ) 


g 4 « 


(49) 


and K is a purely periodic series. 

Thus, we reduce the problem to integration of the partial differential equations: 


3w NI an „ h ? an „ aw 


aw aw aw Wfl M Vll -- U Wll jr ^ If 

?1 T- + g 2 s — + g 3 ^ — + 64 ^ — =N_ + M — K_ 

Bv t Bv 2 Bv 3 * Bv 4 Bu u 2 Bw 2 Bv 4 


^ "if e o sinw i) ’ ( 19t> ) 


ar ar ar ar ,, . ar p 

g > + e 2 — + g 3 — + g 4 — = (1 - + 83 r k 


a^ T 3v2 - *3 3Vj - & 4 3 V4 


h2 h o - . , . „ ar 

+ — • — R C sin (w„ - v 0 ) - K - — 

h 2 U 2 3v 4 


(46a) 


However, the right sides of these equations must not contain constant terms. 
From Equation 49 it can be seen easily that 


K 0 (e 0 k) 


and that it is a small quantity. Numerically the terms, 



and k — 

Bv. 


are at least of the fourth order in Consequently, the solution of the partial differential Equa- 
tions 19b and 46a by means of iteration represents a fast convergent process. 


Some simplifications are possible in the case of an artificial satellite. If we include the ef- 
fects of the tesseral harmonics, then the argument v 4 is the sidereal time on the zero meridian. 
Designating the angular velocity of rotation of the earth by n', we have 


dv 


X 
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Letting g 4 be the constant part in the right side of the last equation and putting 


we have 



_d_ _ 

dv ^ 1 3v 


3 3 x o v o 

+ ^2 ^ ' + C 6*? - §4 ) ^ — K - 


3v, 


3v, 


because the arguments, v 3 and v 4 appear in the form of the difference v 3 - v 4 . Evidently, 


KcO 



If we consider the case of geodetic satellites with small eccentricities, then the quantity n' e 0 /h 3 0 
is small enough to justify the use of iteration. For example, for the satellite with 


a a 1.87 R and e = 0.01 


the factor, n'e 0 /hg is of the order 10~ 3 . 

If only the zonal harmonics are considered, then, 


d 3 3 

= Bl 3^ + 3^ 


and Equations 19b and 46a become 


aw aw 9n h o an . /aw h o \ 

, _ + g, _ = + M— — + (1 - g 2 ) — e 0 sin w 

1 3vj 2 3 v 2 u 2 ^w 2 \^ w i ^ / 


and 


3p 3P /(| . 3P - r h 2 h o - . , , 

Si + g 2 = + g 3 rx k + — • — KCsin (W 2 - v 2 ). 


3vj 2 3 v 2 


K u 2 


CONSTANTS OF INTEGRATION 


Determining integration constants in this theory is simpler than in the author's previous 
theory (Reference 5). 
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In the lunar case, the series for k 11 , k 22 , and ^ 32 are the cosine series and, consequently, 
they contain the additive constant of integration. The series for ^ J2 , k 2l , and k 31 are the sine 
series and they do not contain any such constants. The same conclusions are valid also for arti- 
ficial satellites. We conclude that the form of the additive constant of integration in r i is, 


(1 + c[) cos w 2 ; 


in r 2 it is 


(cos i Q + c') sin w 2 t 


and in r 3 it is 


(sin i Q + c') sin w 2 - 

We shall determine c 2 and c 3 in such a way that cos i 0 sin v 2 will be the only term in P 2> which is 
of the zero order, and sin i 0 sin v 2 will be the only term of zero order in r 3 . The series for r t 
will be the cosine series in v n v 2 , v 3 , and v 4 , and P 2 and P 3 will be the sine series in these argu- 
ments. We conclude that the term in 


+ (1 - g 2 ) — + g r 
2 Bw 2 + g3 2 


independent from v 1? v 2 , v 3 , and v 4 must contain sin w 2 as a factor. In a similar way, the terms in 


c>r 2 


and in 


3r 3 


which are independent from v 1? v 2 , v 3 , and v 4 must have cosw 2 as a factor. In order to avoid the 
secular terms in the components of r, we must remove the terms containing the argument w 2 alone 
from the derivatives of r. Letting be the coefficient of sin w 2 in 




sin (w, - v, - ), 
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K 2 be the coefficient of cos w 2 in 


h 2 h o r . . , . 

-T C ^23 Sin < W 2 - V 2>’ 

h 0 U 


and be the coefficient of cosw 2 in 


we have to set 


h 2 h o - . . , . 

— * — C s m (w 0 - v 9 ) 

h 2 u 2 33 


- (1 - g 2 ) (1 + cj) + g 3 (cos i 0 + cj) + Kj = 0 


(1 - g 2 ) (cos i 0 + cj) - g 3 (1 + cj) + K 2 = 0 


and 


(50) 

(51) 


(1 - g 2 ) (sin i 0 + cj) + K 3 ~ 0, (52) 

and because \ is of the order of perturbation, is of higher order relative to k 2 and K 3 . 

We determine 1 - g 2 and g 3 from Equations 50 and 51, and Equation 52 serves as check. 

The constants cj , cj , and cj in Equations 50 through 52 can be taken from the previous ap- 
proximation. It remains for us to discuss the determination of cj , cj, and c j. Let [r i ] 9 [r 2 ] , [r 3 ] 
be the values obtained by the formal process of integration, without adding the constants. Then, 


Tj = (1 + cj) cos w 2 + [Pj], 

P 2 = (cos i Q + cj) sin w 2 + [P 2 ] , 


and 


P 3 - (sin i 0 + cj) sin w 2 + [P 3 ] . 


Let IXjJq, [r 2 ] 0 , and [P 3 ] 0 be the values of [rj , [r 2 ], and [P 3 ] obtained for w 2 = 0. 
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The quantity f is a unit vector and it must be 


(1 + c[) 2 + const, part in [[Pj ]g + + [P3] = 1 

From this last equation 1 + C J is obtained without any ambiguity. Determining constants in the 
given theory is simpler and easier than in the theory based on use of \ v \ 2 , \ 3> and . 

Equation 11 can be written in the standard form 


W = 3 + Y cos W A + ¥ sin w x , 


(lla) 


where in our case 



(53) 


r = 


h , n ° 

e - — cos cp - — 
h 0 h 


(54) 


and 


¥ = e — s in <f>. 
"0 


If [3J , [T], [¥] , and [W] are the values obtained by formal integration, then 


(55) 


3 = [3] 


+ C 0 


r = [Y] + cj, 

'V = [V] , 


[W] = [H] + [Y] cos Wj + [V] sin w lf 


and 


W _ [W] + c" + c" cos W r 

where c" and c " are the constants of integration. 
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We have 


[g] = the part of [w] which does not contain 

rri = ( [w] _ [g]} 

w x = 0 


and 


m = { [ W] - CH] > _ 

w i 


From Equation 53 we obtain: 



h 2 



+ * • * * 



(53a) 


and 



(53b) 


h~ n 



+ • • • » 


where g 2 , g 3 ? ... can be taken from the previous approximation. We determine c" in such a way 
that h 2 /h 2 - 1 does not contain a constant term. Then 


c" = - the constant term in 





where c" is of the higher order compared to g* From Equation 12 we obtain 


u 

h 2 





w, 


or 


u 

h 2 


0 u r in i " 11 

= — •— + LWJ + c + c cos v, 
h 2 h 2 


+ 



W 


(56) 


(57) 
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In the last equation (h/h 0 ) _l and W can be taken from the previous approximation. We choose c" in 
such a way that no term of the form A cos v 1 is present except for the term, e 0 cos v,. 

DETERMINATION OF TIME 

We introduce the disturbed time, z, by means of the equation 


1 dv 1 

dl ' 


Introducing the eccentric anomaly, e, by means of the equation, 


r = u -1 = h 0 2 (1 + e 0 cos VjT 1 = a 0 (1 - e 0 cos e) ■ 


(58) 


where 


= h o 2 


(1 




we can write Equation 58 as a Kepler's equation of the form, 


e - e Q sin e r { n Q gj z 


where 


Setting 

Z - 


we obtain from Equations 3, 12a and 58 


dn 0 8 j 


dv 


= (1 - e 


2 \ 3/2 
(V 


-3/2 

a o 


t + hz 


f-Y 

1 - A (1 + ,) 2 

\u / 

L h o J 


Making use of Equation 53b we then obtain, 


dn Q 8z 
dv 


-a-e^(|y [(2 v.|h) 

+ H 2 + 3 V + V 2 ) + . . . .J 
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CONCLUSIONS 


The lunar theory presented here is a numerical one, based on the application of the process 
of iteration. The number of cycles for the satellites of outer planets will be rather small, be- 
cause the final accuracy of the computation for such satellites can be about 10~ 5 in 1 + v and 

0.001° in the angles. The decision concerning the choice of terms in the development is left to the 
machine and thus there is no danger that by accident some influential terms will be omitted. We 
assume that the eccentricity is small or moderate. Then validity of the theory depends upon the 
value of the parallactic factor and upon the presence of the resonance effects. 

Recently Mr. M. Charnow has programmed the author T s version of Hansen's lunar theory 
(Reference 5). He computed for comparison and check the coordinates of the moon using Hansen 
input data. The comparison is turned to be quite satisfactory. Thus, the basic ideas expressed in 
this as well as in the author's previous works are shown to be applicable in practice and we are 
planning to use them to develop the theories of some satellites of outer planets. 
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Appendix A 


N' 


K' 


W 

W 


c o 

1 "Si 

g 2 

g 3 

h o 

1/h 



r 

— » o 

r 

— o 

r 

u 
u * 
u 

V 

V 1 = ^1 v - "o 

v 2 = g 2 V + o ; 0 


List of Symbols 

- periodic part in -1/2 (6 + &) , 

- periodic part in + 1/2 (<9 - &) , 

- Hansen T s function determining the perturbations in the orbital plane, 



- osculating eccentricity of the satellite, 

- mean value of the eccentricity of the satellite, 

- mean motion of x with respect to v, 

- mean motion of the argument of the latitude, 

- mean motion of the ascending node, 

- mean value of h, 

- areal velocity of the satellite, 

- osculating inclination of the orbital plane of the satellite toward the orbital plane 
of the sun, 

- mass of the sun. The orbit of the sun is taken to be elliptic, 

= I r I . 

= in 

- position vector of the satellite, 

- position vector of the sun, 

- unit vector of r, 

- the unit vector of r', 

= 1/r, 

= u' = 1/r' , 

- the mean value of u , 

- true orbital longitude of the satellite, 

- mean true anomaly of the satellite, 

- mean argument of the latitude of the satellite, 
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v 3 = g 3 v -« e Q - mean longitude of the ascending node of the satellite, 

W - true orbital longitude of the fictitious satellite, 

Wj = gi w -? 7 0 - mean true anomaly of the fictitious satellite, 
w 2 - g 2 w + cj 0 - mean argument of the latitude of the fictitious satellite, 
z - pseudo-time (the disturbed time), 

8z = z ~ t - perturbations of time, 

r - unit yector of the fictitious satellite with respect to the system of coordinates as- 
sociated with the mean orbital plane. This plane is defined here as a plane having 
the real inclination i but whose longitude of the ascending node is equal to v 3 , 

r - unit vector of the real satellite with respect to the mean orbital plane, 

n - disturbing function of the satellite. The mass of the planet and the gravitation 
constant are chosen to be one. The mass of the satellite is supposed to be 
negligible, 

n - disturbing function in which the elliptic and the non-elliptic v 2 are separated by 
replacing in n the vector r by the vector r , 

X - true orbital longitude of the pericenter of the satellite, 

0() - mean value of x? 

<t> - periodic part of x> 

0 - longitude of the osculating ascending node, 
a - distance of the departure point from the ascending node, 
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